The prompt and base code for this project can be found on Datacamp. The basis of this is to explore the volatility of Zero Coupon Bond Yields (from year maturities 1-30) and see how that has changed since the 1960’s. Below, I loaded the necessary librares and did some data exploration of the time series data obtained from Quandl. This was obtained making an API call.

Necessary packages

defaultW <- getOption("warn") 

options(warn = -1) 

library(dygraphs)
library(Quandl)
library(dplyr)
library(rugarch)
library(plotly)
library(RColorBrewer)

options(warn = defaultW)
#get dataset, filter for coronavirus dates
df <- Quandl("FED/SVENY", api_key="15cxBvvbCzucYDswsDfJ")
covid_dates <- subset(df, Date > '2020-01-01' & Date < '2020-08-20')

str(df)
## 'data.frame':    14755 obs. of  31 variables:
##  $ Date   : Date, format: "2020-08-14" "2020-08-13" ...
##  $ SVENY01: num  0.157 0.162 0.162 0.159 0.146 ...
##  $ SVENY02: num  0.157 0.171 0.163 0.155 0.135 ...
##  $ SVENY03: num  0.197 0.215 0.202 0.187 0.161 ...
##  $ SVENY04: num  0.26 0.278 0.261 0.24 0.209 ...
##  $ SVENY05: num  0.334 0.351 0.331 0.304 0.268 ...
##  $ SVENY06: num  0.412 0.426 0.404 0.373 0.332 ...
##  $ SVENY07: num  0.491 0.502 0.478 0.443 0.398 ...
##  $ SVENY08: num  0.568 0.577 0.551 0.512 0.463 ...
##  $ SVENY09: num  0.643 0.649 0.621 0.58 0.527 ...
##  $ SVENY10: num  0.715 0.717 0.688 0.644 0.589 ...
##  $ SVENY11: num  0.784 0.783 0.752 0.706 0.649 ...
##  $ SVENY12: num  0.849 0.846 0.813 0.766 0.706 ...
##  $ SVENY13: num  0.911 0.906 0.872 0.823 0.762 ...
##  $ SVENY14: num  0.971 0.964 0.927 0.878 0.815 ...
##  $ SVENY15: num  1.028 1.019 0.98 0.93 0.866 ...
##  $ SVENY16: num  1.082 1.071 1.031 0.98 0.915 ...
##  $ SVENY17: num  1.134 1.121 1.08 1.028 0.962 ...
##  $ SVENY18: num  1.18 1.17 1.13 1.07 1.01 ...
##  $ SVENY19: num  1.23 1.22 1.17 1.12 1.05 ...
##  $ SVENY20: num  1.28 1.26 1.21 1.16 1.09 ...
##  $ SVENY21: num  1.32 1.3 1.25 1.2 1.13 ...
##  $ SVENY22: num  1.36 1.34 1.29 1.24 1.17 ...
##  $ SVENY23: num  1.4 1.38 1.33 1.28 1.21 ...
##  $ SVENY24: num  1.44 1.42 1.37 1.31 1.25 ...
##  $ SVENY25: num  1.48 1.45 1.4 1.35 1.28 ...
##  $ SVENY26: num  1.52 1.49 1.43 1.38 1.31 ...
##  $ SVENY27: num  1.55 1.52 1.46 1.41 1.34 ...
##  $ SVENY28: num  1.58 1.55 1.5 1.44 1.38 ...
##  $ SVENY29: num  1.61 1.58 1.52 1.47 1.4 ...
##  $ SVENY30: num  1.64 1.61 1.55 1.5 1.43 ...
##  - attr(*, "freq")= chr "daily"

Plotting daily estimates for yall zero coupon yields for 2020

dfasxts <- as.xts(x = df[, -1], order.by = df$Date)
covid_dates <- as.xts(x = covid_dates[, -1], order.by = covid_dates$Date)
dygraph(covid_dates, main = "All Zero Coupon Yields (1-30) 2020", ylab = "Value") %>%
            dyAxis('x', axisLabelFontSize = 12) %>%
            dyRangeSelector()
df$Date <- as.Date(df$Date)
df$year <- format(as.Date(df$Date, format="%m/%d/%Y"),"%Y")
df <- select(df, -Date)
df <- na.omit(df)

df <- df %>%
  group_by(year) %>%
  summarise_all(mean)
SVENY01 <- df$SVENY01
SVENY10 <- df$SVENY10
SVENY30 <- df$SVENY30


data <- data.frame(df, SVENY01, SVENY10, SVENY30)
data <- select(data, year, SVENY01, SVENY10, SVENY30)
fig <- plot_ly(data, x = data$year, y = ~SVENY01, name = 'SVENY01', type = 'scatter', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY10, name = 'SVENY10', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY30, name = 'SVENY30', mode = 'lines + markers')
fig <- fig %>% layout(title = "",
         xaxis = list(title = "year"),
         yaxis = list (title = "value"))

1, 10, and 30 year maturity bond yield averages across time

fig

### Bond yield heat map

row.names(df) <- df$year
## Warning: Setting row names on a tibble is deprecated.
df <- select(df, -year)
df_matrix <- data.matrix(df)

df_heatmap <- heatmap(df_matrix, Rowv=NA, Colv=NA, col = brewer.pal(6, "Greys"), scale="column", margins=c(4,1), main = "Heatmap")

# plotting the evaluation of bond yields
library(viridisLite)
yields <- dfasxts
plot.type <- "single"
plot.palette <- magma(n = 30)
asset.names <- colnames(dfasxts)
plot.zoo(x = dfasxts, plot.type = "single", col = plot.palette, ylab = "", xlab = "")
legend(x = "topleft", legend = asset.names, col = plot.palette, cex = 0.45, lwd = 3)

dfasxts_d <- diff(dfasxts)

plot.zoo(x = dfasxts_d, plot.type = "multiple", ylim = c(-0.5, 0.5), cex.axis = 0.7, ylab = 1:30, col = plot.palette, main = "", xlab = "")

dfasxts <- dfasxts_d["2000/",]
x_1 <- dfasxts[,"SVENY01"]
x_20 <- dfasxts[, "SVENY20"]

# Plot the autocorrelations of the yield changes)
par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,2))
acf_1 <- acf(x_1)
acf_20 <- acf(x_20)

# Plot the autocorrelations of the absolute changes of yields
acf_abs_1 <- acf(abs(x_1))
acf_abs_20 <- acf(abs(x_20))

spec <- ugarchspec(distribution.model = "sstd")


fit_1 <- ugarchfit(x_1, spec = spec)


vol_1 <- sigma(fit_1)
res_1 <- scale(residuals(fit_1, standardize = TRUE)) * sd(x_1) + mean(x_1)


merge_1 <- merge.xts(x_1, vol_1, res_1)
plot.zoo(merge_1, xlab = "Year")

####################

fit_20 <- ugarchfit(x_20, spec = spec)


vol_20 <- sigma(fit_20)
res_20 <- scale(residuals(fit_20, standardize = TRUE)) * sd(x_20) + mean(x_20)


merge_20 <- merge.xts(x_20, vol_20, res_20)
plot.zoo(merge_20, xlab = "Year")

par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,1))
hist(res_1)
hist(res_20)

ugarchspec()
## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : sGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(1,0,1)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE
density_x_1 <- density(x_1)
density_res_1 <- density(res_1)


plot(density_x_1)
lines(density_res_1, col = "red")


norm_dist <- dnorm(seq(-0.4, 0.4, by = .01), mean = mean(x_1), sd = sd(x_1))
lines(seq(-0.4, 0.4, by = .01), 
      norm_dist, 
      col = "darkblue"
     )

# Add legend
legend <- c("Before GARCH", "After GARCH", "Normal distribution")
legend("topleft", legend = legend, 
       col = c("black", "red", "darkblue"), lty=c(1,1))

distribution <- qnorm


qqnorm(x_1, ylim = c(-0.5, 0.5))
qqline(x_1, distribution = distribution, col = "darkgreen")


par(new=TRUE)
qqnorm(res_1 * 0.614256270265139, col = "red", ylim = c(-0.5, 0.5))
qqline(res_1 * 0.614256270265139, distribution = distribution, col = "darkgreen")
legend("topleft", c("Before GARCH", "After GARCH"), col = c("black", "red"), pch=c(1,1))

Conclusion

Garch modeling brought the residuals closer to normal distribution. Year 1 yield deviates more from a normally distributed white noise process than a 20 year yield.